Homogeneous nucleation rate of methane hydrate formation under experimental conditions from seeding simulations.

In this work, we shall estimate via computer simulations the homogeneous nucleation rate for the methane hydrate at 400 bars for a supercooling of about 35 K. The TIP4P/ICE model and a Lennard-Jones center were used for water and methane, respectively. To estimate the nucleation rate, the seeding technique was employed. Clusters of the methane hydrate of different sizes were inserted into the aqueous phase of a two-phase gas-liquid equilibrium system at 260 K and 400 bars. Using these systems, we determined the size at which the cluster of the hydrate is critical (i.e., it has 50% probability of either growing or melting). Since nucleation rates estimated from the seeding technique are sensitive to the choice of the order parameter used to determine the size of the cluster of the solid, we considered several possibilities. We performed brute force simulations of an aqueous solution of methane in water in which the concentration of methane was several times higher than the equilibrium concentration (i.e., the solution was supersaturated). From brute force runs, we infer the value of the nucleation rate for this system rigorously. Subsequently, seeding runs were carried out for this system, and it was found that only two of the considered order parameters were able to reproduce the value of the nucleation rate obtained from brute force simulations. By using these two order parameters, we estimated the nucleation rate under experimental conditions (400 bars and 260 K) to be of the order of log10 (J/(m3 s)) = -7(5).

Maximum in density of electrolyte solutions: Learning about ion-water interactions and testing the Madrid-2019 force field.

In this work, we studied the effect of Li+, Na+, K+, Mg2+, and Ca2+ chlorides and sulfates on the temperature of maximum density (TMD) of aqueous solutions at room pressure. Experiments at 1 molal salt concentration were carried out to determine the TMD of these solutions. We also performed molecular dynamics simulations to estimate the TMD at 1 and 2 m with the Madrid-2019 force field, which uses the TIP4P/2005 water model and scaled charges for the ions, finding an excellent agreement between experiment and simulation. All the salts studied in this work shift the TMD of the solution to lower temperatures and flatten the density vs temperature curves (when compared to pure water) with increasing salt concentration. The shift in the TMD depends strongly on the nature of the electrolyte. In order to explore this dependence, we have evaluated the contribution of each ion to the shift in the TMD concluding that Na+, Ca2+, and SO4 2- seem to induce the largest changes among the studied ions. The volume of the system has been analyzed for salts with the same anion and different cations. These curves provide insight into the effect of different ions upon the structure of water. We claim that the TMD of electrolyte solutions entails interesting physics regarding ion-water and water-water interactions and should, therefore, be considered as a test property when developing force fields for electrolytes. This matter has been rather unnoticed for almost a century now and we believe it is time to revisit it.

Phase diagram of the NaCl-water system from computer simulations.

NaCl aqueous solutions are ubiquitous. They can crystallize into ice, NaCl, or NaCl · 2H2O depending on the temperature-concentration conditions. These crystallization transitions have important implications in geology, cryopreservation, or atmospheric science. Computer simulations can help understand the crystallization of these solids, which requires a detailed knowledge of the equilibrium phase diagram. We use molecular simulations in which we put at contact the solution with the solid of interest to determine points of the solid-solution coexistence lines. We follow two different approaches, one in which we narrow down the melting temperature for a given concentration and the other in which we equilibrate the concentration for a given temperature, obtaining consistent results. The phase diagram thus calculated for the selected model (TIP4P/2005 for water molecules and Joung-Cheatham for the ions) correctly predicts coexistence between the solution and ice. We were only able to determine NaCl · 2H2O-solution coexistence points at higher temperatures and concentrations than in the experiment, so we could not establish a direct comparison in this case. On the other hand, the model underestimates the concentration of the solution in equilibrium with the NaCl solid. Our results, alongside other literature evidence, seem to indicate that ion-ion interactions are too strong in the model. Our work is a good starting point for the improvement of the potential model and for the study of the nucleation kinetics of the solid phases involved in the phase diagram.

The Young-Laplace equation for a solid-liquid interface.

The application of the Young-Laplace equation to a solid-liquid interface is considered. Computer simulations show that the pressure inside a solid cluster of hard spheres is smaller than the external pressure of the liquid (both for small and large clusters). This would suggest a negative value for the interfacial free energy. We show that in a Gibbsian description of the thermodynamics of a curved solid-liquid interface in equilibrium, the choice of the thermodynamic (rather than mechanical) pressure is required, as suggested by Tolman for the liquid-gas scenario. With this definition, the interfacial free energy is positive, and the values obtained are in excellent agreement with previous results from nucleation studies. Although, for a curved fluid-fluid interface, there is no distinction between mechanical and thermal pressures (for a sufficiently large inner phase), in the solid-liquid interface, they do not coincide, as hypothesized by Gibbs.

Temperature of maximum density and excess thermodynamics of aqueous mixtures of methanol.

In this work, we present a study of representative excess thermodynamic properties of aqueous mixtures of methanol over the complete concentration range, based on extensive computer simulation calculations. In addition to test various existing united atom model potentials, we have developed a new force-field which accurately reproduces the excess thermodynamics of this system. Moreover, we have paid particular attention to the behavior of the temperature of maximum density (TMD) in dilute methanol mixtures. The presence of a temperature of maximum density is one of the essential anomalies exhibited by water. This anomalous behavior is modified in a non-monotonous fashion by the presence of fully miscible solutes that partly disrupt the hydrogen bond network of water, such as methanol (and other short chain alcohols). In order to obtain a better insight into the phenomenology of the changes in the TMD of water induced by small amounts of methanol, we have performed a new series of experimental measurements and computer simulations using various force fields. We observe that none of the force-fields tested capture the non-monotonous concentration dependence of the TMD for highly diluted methanol solutions.

Reverse Monte Carlo modeling in confined systems.

An extension of the well established Reverse Monte Carlo (RMC) method for modeling systems under close confinement has been developed. The method overcomes limitations induced by close confinement in systems such as fluids adsorbed in microporous materials. As a test of the method, we investigate a model system of (36)Ar adsorbed into two zeolites with significantly different pore sizes: Silicalite-I (a pure silica form of ZSM-5 zeolite, characterized by relatively narrow channels forming a 3D network) at partial and full loadings and siliceous Faujasite (which exhibits relatively wide channels and large cavities). The model systems are simulated using grand canonical Monte Carlo and, in each case, its structure factor is used as input for the proposed method, which shows a rapid convergence and yields an adsorbate microscopic structure in good agreement with that of the model system, even to the level of three body correlations, when these are induced by the confining media. The application to experimental systems is straightforward incorporating factors such as the experimental resolution and appropriate q-sampling, along the lines of previous experiences of RMC modeling of powder diffraction data including Bragg and diffuse scattering.

Free energy calculations for molecular solids using GROMACS.

In this work, we describe a procedure to evaluate the free energy of molecular solids with the GROMACS molecular dynamics package. The free energy is calculated using the Einstein molecule method that can be regarded as a small modification of the Einstein crystal method. Here, the position and orientation of the molecules is fixed by using an Einstein field that binds with harmonic springs at least three non-collinear atoms (or points of the molecule) to their reference positions. The validity of the Einstein field is tested by performing free-energy calculations of methanol, water (ice), and patchy colloids molecular solids. The free energies calculated with GROMACS show a very good agreement with those obtained using Monte Carlo and with previously published results.

Three-dimensional patchy lattice model for empty fluids.

The phase diagram of a simple model with two patches of type A and ten patches of type B (2A10B) on the face centred cubic lattice has been calculated by simulations and theory. Assuming that there is no interaction between the B patches the behavior of the system can be described in terms of the ratio of the AB and AA interactions, r. Our results show that, similarly to what happens for related off-lattice and two-dimensional lattice models, the liquid-vapor phase equilibria exhibit reentrant behavior for some values of the interaction parameters. However, for the model studied here the liquid-vapor phase equilibria occur for values of r lower than 1/3, a threshold value which was previously thought to be universal for 2AnB models. In addition, the theory predicts that below r=1/3 (and above a new condensation threshold which is <1/3) the reentrant liquid-vapor equilibria are so extreme that it exhibits a closed loop with a lower critical point, a very unusual behavior in single-component systems. An order-disorder transition is also observed at higher densities than the liquid-vapor equilibria, which shows that the liquid-vapor reentrancy occurs in an equilibrium region of the phase diagram. These findings may have implications in the understanding of the condensation of dipolar hard spheres given the analogy between that system and the 2AnB models considered here.

Heat capacity of water: A signature of nuclear quantum effects.

In this note we present results for the heat capacity at constant pressure for the TIP4PQ/2005 model, as obtained from path-integral simulations. The model does a rather good job of describing both the heat capacity of ice I(h) and of liquid water. Classical simulations using the TIP4P/2005, TIP3P, TIP4P, TIP4P-Ew, simple point charge/extended, and TIP5P models are unable to reproduce the heat capacity of water. Given that classical simulations do not satisfy the third law of thermodynamics, one would expect such a failure at low temperatures. However, it seems that for water, nuclear quantum effects influence the heat capacities all the way up to room temperature. The failure of classical simulations to reproduce C(p) points to the necessity of incorporating nuclear quantum effects to describe this property accurately.

Quantum effects on the maximum in density of water as described by the TIP4PQ/2005 model.

Path integral simulations have been performed to determine the temperature of the maximum in density of water of the rigid, nonpolarizable TIP4PQ/2005 model treating long range Coulombic forces with the reaction field method. A maximum in density is found at 280 K, just 3 K above the experimental value. In tritiated water the maximum occurs at a temperature about 12 K higher than in water, in reasonable agreement with the experimental result. Contrary to the usual assumption that the maximum in classical water is about 14 K above that in water, we found that for TIP4PQ/2005 this maximum is about 30 K above. For rigid water models the internal energy and the temperature of maximum density do not follow a linear behavior when plotted as a function of the inverse of the hydrogen mass. In addition, it is shown that, when used with Ewald sums, the TIP4PQ/2005 reproduces quite nicely not only the maximum in density of water, but also the liquid densities, the structure of liquid water and the vaporization enthalpy. It was shown in a previous work that it also reproduces reasonably well the density and relative stabilities of ices. Therefore the TIP4PQ/2005 model, while still simple, allows one to analyze the interplay between quantum effects related to atomic masses and intermolecular forces in water.

The phase diagram of water at high pressures as obtained by computer simulations of the TIP4P/2005 model: the appearance of a plastic crystal phase.

In this work the high pressure region of the phase diagram of water has been studied by computer simulation by using the TIP4P/2005 model of water. Free energy calculations were performed for ices VII and VIII and for the fluid phase to determine the melting curve of these ices. In addition, molecular dynamics simulations were performed at high temperatures (440 K) observing the spontaneous freezing of the liquid into a solid phase at pressures of about 80,000 bar. The analysis of the structure obtained lead to the conclusion that a plastic crystal phase was formed. In the plastic crystal phase the oxygen atoms were arranged forming a body center cubic structure, as in ice VII, but the water molecules were able to rotate almost freely. Free energy calculations were performed for this new phase, and it was found that for TIP4P/2005 this plastic crystal phase is thermodynamically stable with respect to ices VII and VIII for temperatures higher than about 400 K, although the precise value depends on the pressure. By using Gibbs-Duhem simulations, all coexistence lines were determined, and the phase diagram of the TIP4P/2005 model was obtained, including ices VIII and VII and the new plastic crystal phase. The TIP4P/2005 model is able to describe qualitatively the phase diagram of water. It would be of interest to study if such a plastic crystal phase does indeed exist for real water. The nearly spherical shape of water makes possible the formation of a plastic crystal phase at high temperatures. The formation of a plastic crystal phase at high temperatures (with a bcc arrangements of oxygen atoms) is fast from a kinetic point of view occurring in about 2 ns. This is in contrast to the nucleation of ice Ih which requires simulations of the order of hundreds of ns.

Computing the free energy of molecular solids by the Einstein molecule approach: ices XIII and XIV, hard-dumbbells and a patchy model of proteins.

The recently proposed Einstein molecule approach is extended to compute the free energy of molecular solids. This method is a variant of the Einstein crystal method of Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. In order to show its applicability, we have computed the free energy of a hard-dumbbell solid, of two recently discovered solid phases of water, namely, ice XIII and ice XIV, where the interactions between water molecules are described by the rigid nonpolarizable TIP4P/2005 model potential, and of several solid phases that are thermodynamically stable for an anisotropic patchy model with octahedral symmetry which mimics proteins. Our calculations show that both the Einstein crystal method and the Einstein molecule approach yield the same results within statistical uncertainty. In addition, we have studied in detail some subtle issues concerning the calculation of the free energy of molecular solids. First, for solids with noncubic symmetry, we have studied the effect of the shape of the simulation box on the free energy. Our results show that the equilibrium shape of the simulation box must be used to compute the free energy in order to avoid the appearance of artificial stress in the system that will result in an increase in the free energy. In complex solids, such as the solid phases of water, another difficulty is related to the choice of the reference structure. As in some cases there is no obvious orientation of the molecules; it is not clear how to generate the reference structure. Our results will show that, as long as the structure is not too far from the equilibrium structure, the calculated free energy is invariant to the reference structure used in the free energy calculations. Finally, the strong size dependence of the free energy of solids is also studied.

Properties of ices at 0 K: a test of water models.

The properties of ices Ih, II, III, V, and VI at zero temperature and pressure are determined by computer simulation for several rigid water models (SPC/E, TIP5P, TIP4P/Ice, and TIP4P/2005). The energies of the different ices at zero temperature and pressure (relative to the ice II energy) are compared to the experimental results of Whalley [J. Chem. Phys. 81, 4087 (1984)]. TIP4P/Ice and TIP4P/2005 provide a qualitatively correct description of the relative energies of the ices at these conditions. In fact, only these two models provide the correct ordering in energies. For the SPC/E and TIP5P models, ice II is the most stable phase at zero temperature and pressure whereas for TIP4P/Ice and TIP4P/2005 ice Ih is the most stable polymorph. These results are in agreement with the relative stabilities found at higher temperatures. The solid-solid phase transitions at 0 K are determined. The predicted pressures are in good agreement with those obtained from free energy calculations.

Phase diagram of model anisotropic particles with octahedral symmetry.

The phase diagram for a system of model anisotropic particles with six attractive patches in an octahedral arrangement has been computed. This model for a relatively narrow value of the patch width where the lowest-energy configuration of the system is a simple cubic crystal. At this value of the patch width, there is no stable vapor-liquid phase separation, and there are three other crystalline phases in addition to the simple cubic crystal that is most stable at low pressure. First, at moderate pressures, it is more favorable to form a body-centered-cubic crystal, which can be viewed as two interpenetrating, and almost noninteracting, simple cubic lattices. Second, at high pressures and low temperatures, an orientationally ordered face-centered-cubic structure becomes favorable. Finally, at high temperatures a face-centered-cubic plastic crystal is the most stable solid phase.

A density-functional study of the structures, binding energies and total spins of Ni-Fe clusters using nonlocal norm-conserving pseudopotentials and the generalized gradient approximation.

We report ab initio calculations of the structures, binding energies, and total spins of the clusters Ni(13), Ni(19), Ni(23), Ni(26), Ni(12)Fe, Ni(11)Fe(2), Ni(18)Fe, Ni(17)Fe(2), Ni(22)Fe, Ni(20)Fe(3), and Ni(25)Fe using a density-functional method that employs linear combination of atomic orbitals as basis sets, nonlocal norm-conserving pseudopotentials, and the generalized gradient approximation for exchange and correlation. Our results show that the Fe-doped Ni clusters, which have icosahedral or polyicosahedral ground-state structures similar to those of the corresponding pure Ni clusters, are most stable with the Fe atoms occupying internal positions, as has also been inferred from experimental results on the adsorption of molecular nitrogen on the cluster surfaces. We also rule out the possibility that the experimentally observed difference between the (nonpolyicosahedral) configurations of N(2)-saturated Ni(26) and N(2)-saturated Ni(25)Fe be due to the influence of the Fe atom on the energy of the underlying metal cluster.

Geometric structure and electronic properties of neutral and anionic Fe2C3 and Fe2C4 clusters, as obtained by density-functional calculations.

We calculated the geometrical structures and electronic properties of neutral and anionic Fe2Cn clusters (n = 3,4) using a density-functional method that employs linear combinations of atomic orbitals as basis sets, standard nonlocal norm-conserving pseudopotentials, and the generalized gradient approximation to exchange and correlation. We show that the ground-state structures of Fe2C3 and Fe2C4 are essentially the same in the neutral and anionic states, namely, planar rings that feature nonadjacent Fe atoms. For the anionic clusters, these findings contrast with previously published results.